Theory of electric polarization induced by inhomogeneity in crystals 
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We develop a general theory of electric polarization induced by inhomogeneity in crystals. We 
show that contributions to polarization can be classified in powers of the gradient of the order 
parameter. The zeroth order contribution reduces to the well-known result obtained by King- 
Smith and Vanderbilt for uniform systems. The first order contribution, when expressed in a two- 
point formula, takes the Chern-Simons 3-form of the vector potentials derived from the Bloch wave 
functions. Using the relation between polarization and charge density, we demonstrate our formula 
by studying charge fractionalization in a two-dimensional dimer model recently proposed. 
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Electric polarization is a fundamental quantity in con- 
densed matter physics, essential to any proper descrip- 
tion of dielectric phenomena of matter. Theoretically, it 
is well established that only the change in polarization 
has physical meaning and it can be quantified by using 
the Berry phase of the electronic wave functions [l|, |j, Isj . 
In practice, the Berry-phase formula is usually expressed 
in terms of the Bloch orbitals. It has been very successful 
in first-principles studies of dielectric properties of oxides 
and other insulating materials. 

While the existing formulation is adequate in periodic 
insulators, a theory of polarization for inhomogeneous 
crystals would find numerous important applications; for 
example, in a class of recently discovered multiferroics, 
the appearance of electric polarization is always accom- 
panied by long- wavelength magnetic structures [J, S l6| . 
A number of phenomenological and microscopic theories 
have been proposed to understand this magnetically in- 
duced polarization [y, 0, |^, M, \1^ ', however, quantitative 
studies of this type of problem still remain in a primitive 
state. The fundamental difficulty lies in the fact that the 
inhomogeneous ordering breaks the translational symme- 
try of the crystal so that Bloch's theorem does not apply. 

In this Letter we present a general framework to calcu- 
late electric polarization in crystals with inhomogeneous 
ordering. Our theory is based on the elementary relation 
between the change in polarization and integrated bulk 
current [2, 13(1 . The latter can be evaluated using the semi- 
classical formalism of Bloch electron dynamics [ll|. We 
find that, in addition to the contribution previously ob- 
tained for uniform systems |1|] , the polarization contains 
an extra contribution proportional to the gradient of the 
order parameter. This extra contribution is expressed 
using the second Chern form of the Berry curvatures de- 
rived from the local Bloch functions. It can also be recast 
into a two-point formula, which depends only on the ini- 
tial and final states, up to an uncertainty quantum after 
spatial averaging. We identify this quantum as the sec- 
ond Chern number in appropriate units. In addition. 



several general conditions for the inhomogeneity-induced 
polarization to be nonzero are also derived. 

To demonstrate our theory, we apply our formula to 
study the problem of charge fractionalization in a two- 
dimensional dimer model recently proposed 13, Il3l | . We 
show that in this model fractional charge appears as a 
result of the ferroelectric domain walls. By using the re- 
lation between polarization and charge density, we calcu- 
late the total charge carried by a vortex in the dimeriza- 
tion pattern and compare it to previous results [ij, ll3| . 
Our approach has the advantage that it can be easily 
incorporated in a band calculation, while previously one 
relied on spectral analysis of the Dirac Hamiltonian per- 
formed in the continuum limit 1^, [l3| . 

General formulation. — Suppose we have an insulating 
crystal with an order parameter m{r) that varies slowly 
in space. We assume that, at least on the mean-field level, 
m{r) can be treated as an external field that couples to 
an operator in the Hamiltonian Ti. Thus, we can formally 
write Ti[m{r)]. As was emphasized in previous work [l|, 
|2, l3|, only the change in polarization P between two 
different states has meaning, and it is given by [14] 



dtjir,t) 



(1) 



where j{r,t) is the bulk current density as the system 
adiabatically evolves from the initial state (i = 0) to the 
final state {t = T). In other words, we assume that the 
two states are connected through a continuous transfor- 
mation of the Hamiltonian 'H[m{r); A] parameterized by 
a scalar A with A(0) = and A(r) = 1. 

In order to find the current density j{r,t), we adopt 
the formalism of semiclassical dynamics of Bloch elec- 
trons [III, which is a powerful tool to investigate the 



influence of slowly varying perturbations on electron dy- 
namics. Within this approach, each electron is described 
by a narrow wave packet localized around Tc and fee h^ 
the phase space. If m{r) varies smoothly compared to 
the width of the wave packet, it is sufficient to study a 



family of local Hamiltonians 7ic[»Ti(rc); A] which assumes 
a fixed value of the order parameter m{rc) in the vicinity 
of r^. Since 7ic['Ti('"c); A] maintains the periodicity of the 
unperturbed crystal, its eigenstates have the Bloch form: 
|V'n(fe,r-c;A)) = e''""^|u„(fc,rc;A)), where |u„(fe,rc;A)) 
is the cell-periodic part of the Bloch functions. Note 
that the rc-dependence of |u„(fc,rc;A)) enters through 
m(rc). We can then expand the wave packet using these 
local Bloch functions. For simplicity, in the following 
derivation we shall confine ourselves to the case of non- 
degenerate bands and hence omit the band index n. 

It has been previously shown that the wave packet cen- 
ter satisfies the following equations of motion (hereafter 
the subscript c on fee and Tc is dropped) [lli | 
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where e is the electron energy and we have introduced 
the notation V^ = d/dka and VJ^ = d/dra- Summation 
over repeated indices is implied throughout our deriva- 
tion. Here, J7 is the Berry curvature obtained from the 
vector potential A. derived from |u(fe,r. A)). For exam- 
ple. 
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Other Berry curvatures are similarly defined. It is note- 
worthy that although the vector potential A depends on 
the phase choice of the wave function |M(fe,r,A)), the 
Berry curvature Jl is a well-defined gauge-invariant quan- 
tity in the parameter space (fc, r, A). 

We now turn to the derivation of P using Eq. ([T]). The 
electronic contribution to polarization is given by 



P = -e / dk dtD{k,r)r 
Jbz Jo 



(5) 



where — e is the electron charge, and D{r, k) is the elec- 
tron density of states, which is modified from its usual 
value of l/(27r)'' in the presence of the Berry curvature, 
D{k,r) ^il + nt^)/{27rr 15]. 

We can solve Tq, from Eq. ((2]) then insert it into Eq. ([5]). 
Collecting terms proportional to A and keeping those up 
to first order in the gradient, we obtain |16[ 



p = p(0)+p(l) ^ 

where P^"' is the zeroth order contribution 

dk 



p(o) = 



/BZ (27r)- Jo 
and P(^^ is the first order contribution 
dk 
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TABLE I: Comparison between P<°' and P'^' 



p{0) 

p{i) 



Two-point formula 



Uncertain quantum 



Chern- Simons 1-form 
Chern-Simons 3-form 



First Chern number 
Second Chern number 



These are the central results of this work. We note that 
P*^"^ has been obtained by King-Smith and Vanderbilt 
for uniform systems [l|, whereas P^^\ being proportional 
to the gradient of m{r), only exists in inhomogeneous 
crystals. 

Two remarks are in order: firstly, although in the 
above derivation we have assumed an inhomogeneous or- 
der parameter, it is obvious that our theory is also ap- 
plicable when the system is subject to a perturbation of 
a spatially-varying external field; secondly, we have only 
considered the electronic contribution to P here. When 
comparing with experiment, one should also include the 
ionic contribution, which is relatively easy to calculate 
because of its classical nature. 

Two-point formula. — We first show that P^^> has the 
desired property that it depends only on the initial and 
final states. The gauge-invariance of Eq. ([5]) allows us to 
evaluate it with any gauge choice. In order to carry out 
the integration over A, we choose the path-independent 
gauge by requiring that the phase difference between 
|u(fc,r,A)) and \u[k + G,r,X)) does not depend on A, 
where G is a reciprocal lattice vector [3|. Under this 
gauge, Eq. ^ can be recast as [l7| 



p(i) ^ 



BZ 
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. (^) 
We recognize that the integrand in the above equation is 

nothing but the Chern-Simons 3-form. 

However, we have paid a price for performing the A- 

integration; namely, the spatially averaged polarization 



{Pa ) = (l/V) J drPa resulting from this two-point 
formula ^ can only be determined modulo a quantum. 
To find the size of the quantum, we consider a cyclic 
change in A. Let us now assume that the order parame- 
ter m{r) is periodic in r. The integral in Eq. ([8]) (after 
a spatial integration) over a closed manifold spanned by 
{ka,k^ rp, A) is an integer called the second Chern num- 
ber [l8|. Since Eq. ^ does not track the evolution of A, 
there is no information of how many cycles A has gone 
through. This is the reason why {Pa ) using Eq. ([9]) 
can only be determined modulo a quantum. Assuming 
m{r) depends on y, we obtain the quantum for Px in 
a three-dimensional system: 



A(Pii)) = 



tydz 



(10) 



where ly is the period of m{y) and a^ is the lattice con- 
stant along z. 



Similarly, the zeroth order contribution P^^^ can also 
be cast into a two-point formula and the uncertain quan- 
tum is given by e/(ayaz) yj. First-principles calculations 
show that in real materials P^"^ is usually smaller than 
this quantum. Hence the ratio between P'^*'^ and P'^"'^-' is 
roughly on the order of ly/ciy. The similarities between 
P^^^f and P^^^f are summarized in Table ID 

Minimal conditions for a finite P^^' . — We now eval- 
uate Eq. ^ using a particular path of A. We write 
7i[m(r);A] = 7i[Am(r)] so that A acts like a "switch" 
of the order m{r), i.e., when A = the system is order- 
less and when A = 1 the order is fully developed. Using 
the relation V^; = VJ^m^V™ and V^ = (m^/A)V™, we 
can recast Eq. ^ as 
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(11) 



As we shall see below, this equation is very useful in 
assessing the general properties of P^), 

Beside having the crystal be inhomogeneous, there are 
three general conditions for P^) ^ be nonzero accord- 
ing to Eq. (fTTj) : (i) the system must be two-dimensional 
or higher; (ii) the order parameter m{r) must have two 
or more components; and (iii) the wave function must 
depend on four or more independent parameters. These 
conditions can be obtained by realizing that the inte- 
grand in Eq. pT|) is actually the second Chern 4-form 
i7 A 51 given in its local expression with respect to the co- 
ordinates {ka, kf3, m^, m^). It is antisymmetric in ka and 
kp, and in TTi^ and niu^ hence condition (i) and (ii). Con- 
dition (iii) follows from the fact that all 4-forms vanish 
identically in three or less dimensions. Based on con- 
dition (iii) we can further deduce that dim(7i) > 2. If 
dim(7i) = 2, 7i has four components. However, since 
shifting and scaling energy has no effect on wave func- 
tions, the wave function can depend on only two indepen- 
dent parameters (for example, the spherical coordinates 
on a 2-sphere 5^) and P^) vanishes in this case. This 
set of conditions puts powerful constraints on possible 
microscopic models that display finite P^). Conditions 
(i) and (iii) can also be obtained directly from Eq. ([H]). 

Let us consider a two-dimensional "minimal" model 
and assume that both the space of m{r) and coordinate 
space are two-dimensional. Because of its antisymmet- 
ric properties, we can write the integrand of Eq. pip as 
^apf-iiyX- Then Eq. (fTTj) takes the following form 



P^) = ex[(V • m)m - (m • V)m] , 



(12) 



Here Xj as a function of m{r), can be spatial dependent. 
Interestingly, if we identify mir) with the magnetiza- 
tion order parameter Af (r), the above result is consistent 
with the Landau- Ginzburg theory of polarization induced 
by spiral magnetic ordering [7]. However, our result IE 



FIG. 1: (color online). The total charge carried by the ferro- 
electric vortex domain wall m{r)e^ = nix + iniy as a function 
of A/m, where A is the staggered sublattice potential, and 
m is the dimerization order parameter. As A increases, the 
difference between the result from the continuum limit (red 
dashed line) and that based on the band calculation (blue 
solid line) becomes significant. The insert shows the two- 
dimensional dimerized square lattice with vr-flux per plaque- 
tte in the absence of the vortex. The hopping amplitude is 



i(l±i 



along the x, y-direction. 



is a direct consequence of the minimal dimensionality and 
we did not invoke any symmetry analysis. For higher di- 
mensions, one will have to carry out a careful symmetry 
analysis of the magnetic groups of the crystal [19| . 

Degenerate bands. — So far, our derivation is for non- 
degenerate bands. The generalization to degenerate 



2l[. The vector potential 



bands is straightforward [2C 
and Berry curvature become matrix-valued and are de- 
fined by 



(A)mn = (UmKVa|u„) , (13) 

fiab = VaA - VbA - i[Aa,Ab] , (14) 

where a,b G (fc,r,A) and \um) and |u„) are degenerate 
bands. We then need to take the trace of Eqs. ([7]) and 
([5]) for the zeroth and first order contributions to P. The 
two-point formula in Eq. also takes the non-Abelian 
Chern-Simons form. 

Fractional charge. — To demonstrate our theory, we 
consider the problem of charge fractionalization in a re- 
cently proposed two-dimensional dimer model [ij, n3^ . 
shown schematically in the inset of Fig. [TJ Introducing 



7i = CTj (g) az, 74 = 1 (Xi cTj; and 75 = 1 ( 
the Hamiltonian a,s Ti. = hajai where 



ay, we can write 



h — t{coskx 



cos ky, A, m^ 



sin/cr 



2y sin ky) 



(15) 



tA is the staggered sublattice potential, t(l ± rux) and 
t(l ± my) are the dimerized hopping amplitudes along x 
and y direction. We choose the Landau gauge so that 
the effect of the tt flux is represented by alternating signs 



of the hopping amphtudes along adjacent rows. It turns 
out that this model is a minimal one satisfying all our 
three conditions: (i) it is two-dimensional; (ii) the order 
parameter m = (m^jmy) has two components; and (iii) 
h (after scaling) can be mapped onto a unit sphere S"* 
with four independent spherical angles. 

It can be verified that the energy spectrum of this 
Hamiltonian consists of two doubly degenerate levels; 
therefore, the non-Abelian formalism is necessary. The 
Berry curvature has SU(2) symmetry ISJj IS^j l23| ; hence, 
i-**-"' always vanishes since the non-Abelian version of 
Eq. ([7]) has vanishing trace. Thus, we will only consider 
P*^^^ in what follows. 

Suppose there is a vortex in the dimerization pattern: 
namely, rrix + iniy — ra{r)e™^ . According to Eq. p^ 
together with the fact that p{r) = — V • P, this ferroelec- 
tric vortex domain wall will carry a polarization charge 
oi Q = J drp{r) = nm? J^^ d0x [l| , which is in general 
fractional. 

To compare with previous results, we shall first eval- 
uate X in the continuum limit. Expanding the Hamilto- 
nian around the Dirac point (7r/2,7r/2), we find, accord- 
ing to Eq. (fnl), 



X-^n 



dk 



dX- 



AA 



(fc2 



l2A2+A2)5/2 



(16) 



Since at large k the integrand decays as fc ^, we can 
extend the integration range of k to infinity and obtain 



X 



47r 777,2 



fl 



A 



VA2 



(17) 



and the total charge carried by the vortex is given by 



Q^n-{1 



VA2 



(18) 



where n is the winding number. This result agrees with 
the spectral analysis of the Dirac Hamiltonian IJ, |l3| . 

The above derivation provides a simple picture of 
charge fractionalization in this type of system: it is a 
direct consequence of the ferroelectric domain wall, and 
the breaking of the sublattice symmetry (A) allows it 
to be irrational. A detailed report including both ID 
and 2D cases will be reported elsewhere. We also calcu- 
late the total charge based on a band calculation using 
Eq. p5|) . shown in Fig. [TJ As A increases, the devia- 
tion between the band calculation and continuum limit 
becomes significant. 

In summary, we have developed a general theory of 
polarization induced by inhomogeneity in crystals. Our 
result lays the foundation for quantitative studies of this 
type of problem. In connection to multiferroics, the min- 
imal conditions for a finite P*^^^ point to general direc- 
tions to aid in the search for microscopic models. In 
addition, we have illustrated our theory by showing that 



the fractional charge in certain models can be understood 
as the polarization charge accompanying ferroelectric do- 
main walls. 
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